GENESIS CGDYN: large-scale coarse-grained MD simulation with dynamic load balancing for heterogeneous biomolecular systems.

Residue-level coarse-grained (CG) molecular dynamics (MD) simulation is widely used to investigate slow biological processes that involve multiple proteins, nucleic acids, and their complexes. Biomolecules in a large simulation system are distributed non-uniformly, limiting computational efficiency with conventional methods. Here, we develop a hierarchical domain decomposition scheme with dynamic load balancing for heterogeneous biomolecular systems to keep computational efficiency even after drastic changes in particle distribution. These schemes are applied to the dynamics of intrinsically disordered protein (IDP) droplets. During the fusion of two droplets, we find that the changes in droplet shape correlate with the mixing of IDP chains. Additionally, we simulate large systems with multiple IDP droplets, achieving simulation sizes comparable to those observed in microscopy. In our MD simulations, we directly observe Ostwald ripening, a phenomenon where small droplets dissolve and their molecules redeposit into larger droplets. These methods have been implemented in CGDYN of the GENESIS software, offering a tool for investigating mesoscopic biological processes using the residue-level CG models.

Acceleration of generalized replica exchange with solute tempering simulations of large biological systems on massively parallel supercomputer.

Generalized replica exchange with solute tempering (gREST) is one of the enhanced sampling algorithms for proteins or other systems with rugged energy landscapes. Unlike the replica-exchange molecular dynamics (REMD) method, solvent temperatures are the same in all replicas, while solute temperatures are different and are exchanged frequently between replicas for exploring various solute structures. Here, we apply the gREST scheme to large biological systems containing over one million atoms using a large number of processors in a supercomputer. First, communication time on a multi-dimensional torus network is reduced by matching each replica to MPI processors optimally. This is applicable not only to gREST but also to other multi-copy algorithms. Second, energy evaluations, which are necessary for the multistate bennet acceptance ratio (MBAR) method for free energy estimations, are performed on-the-fly during the gREST simulations. Using these two advanced schemes, we observed 57.72 ns/day performance in 128-replica gREST calculations with 1.5 million atoms system using 16,384 nodes in Fugaku. These schemes implemented in the latest version of GENESIS software could open new possibilities to answer unresolved questions on large biomolecular complex systems with slow conformational dynamics.

Implementation of residue-level coarse-grained models in GENESIS for large-scale molecular dynamics simulations.

Residue-level coarse-grained (CG) models have become one of the most popular tools in biomolecular simulations in the trade-off between modeling accuracy and computational efficiency. To investigate large-scale biological phenomena in molecular dynamics (MD) simulations with CG models, unified treatments of proteins and nucleic acids, as well as efficient parallel computations, are indispensable. In the GENESIS MD software, we implement several residue-level CG models, covering structure-based and context-based potentials for both well-folded biomolecules and intrinsically disordered regions. An amino acid residue in protein is represented as a single CG particle centered at the Cα atom position, while a nucleotide in RNA or DNA is modeled with three beads. Then, a single CG particle represents around ten heavy atoms in both proteins and nucleic acids. The input data in CG MD simulations are treated as GROMACS-style input files generated from a newly developed toolbox, GENESIS-CG-tool. To optimize the performance in CG MD simulations, we utilize multiple neighbor lists, each of which is attached to a different nonbonded interaction potential in the cell-linked list method. We found that random number generations for Gaussian distributions in the Langevin thermostat are one of the bottlenecks in CG MD simulations. Therefore, we parallelize the computations with message-passing-interface (MPI) to improve the performance on PC clusters or supercomputers. We simulate Herpes simplex virus (HSV) type 2 B-capsid and chromatin models containing more than 1,000 nucleosomes in GENESIS as examples of large-scale biomolecular simulations with residue-level CG models. This framework extends accessible spatial and temporal scales by multi-scale simulations to study biologically relevant phenomena, such as genome-scale chromatin folding or phase-separated membrane-less condensations.

The inherent flexibility of receptor binding domains in SARS-CoV-2 spike protein.

Spike (S) protein is the primary antigenic target for neutralization and vaccine development for the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). It decorates the virus surface and undergoes large motions of its receptor binding domains (RBDs) to enter the host cell. Here, we observe Down, one-Up, one-Open, and two-Up-like structures in enhanced molecular dynamics simulations, and characterize the transition pathways via inter-domain interactions. Transient salt-bridges between RBDA and RBDC and the interaction with glycan at N343B support RBDA motions from Down to one-Up. Reduced interactions between RBDA and RBDB in one-Up induce RBDB motions toward two-Up. The simulations overall agree with cryo-electron microscopy structure distributions and FRET experiments and provide hidden functional structures, namely, intermediates along Down-to-one-Up transition with druggable cryptic pockets as well as one-Open with a maximum exposed RBD. The inherent flexibility of S-protein thus provides essential information for antiviral drug rational design or vaccine development.

Use of multistate Bennett acceptance ratio method for free-energy calculations from enhanced sampling and free-energy perturbation.

Multistate Bennett acceptance ratio (MBAR) works as a method to analyze molecular dynamics (MD) simulation data after the simulations have been finished. It is widely used to estimate free-energy changes between different states and averaged properties at the states of interest. MBAR allows us to treat a wide range of states from those at different temperature/pressure to those with different model parameters. Due to the broad applicability, the MBAR equations are rather difficult to apply for free-energy calculations using different types of MD simulations including enhanced conformational sampling methods and free-energy perturbation. In this review, we first summarize the basic theory of the MBAR equations and categorize the representative usages into the following four: (i) perturbation, (ii) scaling, (iii) accumulation, and (iv) full potential energy. For each, we explain how to prepare input data using MD simulation trajectories for solving the MBAR equations. MBAR is also useful to estimate reliable free-energy differences using MD trajectories based on a semi-empirical quantum mechanics/molecular mechanics (QM/MM) model and ab initio QM/MM energy calculations on the MD snapshots. We also explain how to use the MBAR software in the GENESIS package, which we call mbar_analysis, for the four representative cases. The proposed estimations of free-energy changes and thermodynamic averages are effective and useful for various biomolecular systems.

Structural and energetic analysis of metastable intermediate states in the E1P-E2P transition of Ca2+-ATPase.

Sarcoplasmic reticulum (SR) Ca2+-ATPase transports two Ca2+ ions from the cytoplasm to the SR lumen against a large concentration gradient. X-ray crystallography has revealed the atomic structures of the protein before and after the dissociation of Ca2+, while biochemical studies have suggested the existence of intermediate states in the transition between E1P⋅ADP⋅2Ca2+ and E2P. Here, we explore the pathway and free energy profile of the transition using atomistic molecular dynamics simulations with the mean-force string method and umbrella sampling. The simulations suggest that a series of structural changes accompany the ordered dissociation of ADP, the A-domain rotation, and the rearrangement of the transmembrane (TM) helices. The luminal gate then opens to release Ca2+ ions toward the SR lumen. Intermediate structures on the pathway are stabilized by transient sidechain interactions between the A- and P-domains. Lipid molecules between TM helices play a key role in the stabilization. Free energy profiles of the transition assuming different protonation states suggest rapid exchanges between Ca2+ ions and protons when the Ca2+ ions are released toward the SR lumen.

Optimized Hydrogen Mass Repartitioning Scheme Combined with Accurate Temperature/Pressure Evaluations for Thermodynamic and Kinetic Properties of Biological Systems.

In recent years, molecular dynamics (MD) simulations with larger time steps have been performed with the hydrogen-mass-repartitioning (HMR) scheme, where the mass of each hydrogen atom is increased to reduce high-frequency motion while the mass of a non-hydrogen atom bonded to a hydrogen atom is decreased to keep the total molecular mass unchanged. Here, we optimize the scaling factors in HMR and combine them with previously developed accurate temperature/pressure evaluations. The heterogeneous HMR scaling factors are useful to avoid the structural instability of amino acid residues having a five- or six-membered ring in MD simulations with larger time steps. It also reproduces kinetic properties, namely translational and rotational diffusions, if the HMR scaling factors are applied to only solute molecules. The integration scheme is tested for biological systems that include soluble/membrane proteins and lipid bilayers for about 200 µs MD simulations in total and give consistent results in MD simulations with both a small time step of 2.0 fs and a large, multiple time step integration with time steps of 3.5 fs (for fast motions) and 7.0 fs (for slower motions). We also confirm that the multiple time step integration scheme used in this study provides more accurate energy conservations than the RESPA/C1 and is comparable to the RESPA/C2 in NAMD. In summary, the current integration scheme combining the optimized HMR with accurate temperature/pressure evaluations can provide stable and reliable MD trajectories with a larger time step, which are computationally more than 2-fold efficient compared to the conventional methods.

Elucidation of interactions regulating conformational stability and dynamics of SARS-CoV-2 S-protein.

The ongoing COVID-19 pandemic caused by the new coronavirus, SARS-CoV-2, calls for urgent developments of vaccines and antiviral drugs. The spike protein of SARS-CoV-2 (S-protein), which consists of trimeric polypeptide chains with glycosylated residues on the surface, triggers the virus entry into a host cell. Extensive structural and functional studies on this protein have rapidly advanced our understanding of the S-protein structure at atomic resolutions, although most of these structural studies overlook the effect of glycans attached to the S-protein on the conformational stability and functional motions between the inactive down and active up forms. Here, we performed all-atom molecular dynamics simulations of both down and up forms of a fully glycosylated S-protein in solution as well as targeted molecular dynamics simulations between them to elucidate key interdomain interactions for stabilizing each form and inducing the large-scale conformational transitions. The residue-level interaction analysis of the simulation trajectories detects distinct amino acid residues and N-glycans as determinants on conformational stability of each form. During the conformational transitions between them, interdomain interactions mediated by glycosylated residues are switched to play key roles on the stabilization of another form. Electrostatic interactions, as well as hydrogen bonds between the three receptor binding domains, work as driving forces to initiate the conformational transitions toward the active form. This study sheds light on the mechanisms underlying conformational stability and functional motions of the S-protein, which are relevant for vaccine and antiviral drug developments.

New parallel computing algorithm of molecular dynamics for extremely huge scale biological systems.

In this paper, we address high performance extreme-scale molecular dynamics (MD) algorithm in the GENESIS software to perform cellular-scale molecular dynamics (MD) simulations with more than 100,000 CPU cores. It includes (1) the new algorithm of real-space nonbonded interactions maximizing the performance on ARM CPU architecture, (2) reciprocal-space nonbonded interactions minimizing communicational cost, (3) accurate temperature/pressure evaluations that allows a large time step, and (4) effective parallel file inputs/outputs (I/O) for MD simulations of extremely huge systems. The largest system that contains 1.6 billion atoms was simulated using MD with a performance of 8.30 ns/day on Fugaku supercomputer. It extends the available size and time of MD simulations to answer unresolved questions of biomacromolecules in a living cell.

Group-based evaluation of temperature and pressure for molecular dynamics simulation with a large time step.

Recently, we proposed novel temperature and pressure evaluations in molecular dynamics (MD) simulations to preserve the accuracy up to the third order of a time step, δt [J. Jung, C. Kobayashi, and Y. Sugita, J. Chem. Theory Comput. 15, 84-94 (2019); J. Jung, C. Kobayashi, and Y. Sugita, J. Chem. Phys. 148, 164109 (2018)]. These approaches allow us to extend δt of MD simulations under an isothermal-isobaric condition up to 5 fs with a velocity Verlet integrator. Here, we further improve the isothermal-isobaric MD integration by introducing the group-based evaluations of system temperature and pressure to our previous approach. The group-based scheme increases the accuracy even using inaccurate temperature and pressure evaluations by neglecting the high-frequency vibrational motions of hydrogen atoms. It also improves the overall performance by avoiding iterations in thermostat and barostat updates and by allowing a multiple time step integration such as r-RESPA (reversible reference system propagation algorithm) with our proposed high-precision evaluations of temperature and pressure. Now, the improved integration scheme conserves physical properties of lipid bilayer systems up to δt = 5 fs with velocity Verlet as well as δt = 3.5 fs for fast motions in r-RESPA, respectively.

A singularity-free torsion angle potential for coarse-grained molecular dynamics simulations.

Conventional torsion angle potentials used in molecular dynamics (MD) have a singularity problem when three bonded particles are collinearly aligned. This problem is often encountered in coarse-grained (CG) simulations. Here, we propose a new form of the torsion angle potential, which introduces an angle-dependent modulating function. By carefully tuning the parameters for this modulating function, our method can eliminate the problematic angle-dependent singularity while being combined with existing models. As an example, we optimized the modulating function of the torsion angle potential for popular CG models of biomolecules based on the statistics over experimental structures deposited in the Protein Data Bank. By applying our method to designed and natural biomolecules, we show that the new torsion angle potential is able to eliminate the singularity problem while maintaining the structural features in the original models. Furthermore, by comparing our design with previous methods, we found that our new potential has advantages in computational efficiency and numerical stability. We strongly recommend the usage of our new potential in the CG simulations of flexible molecules.

Scaling molecular dynamics beyond 100,000 processor cores for large-scale biophysical simulations.

The growing interest in the complexity of biological interactions is continuously driving the need to increase system size in biophysical simulations, requiring not only powerful and advanced hardware but adaptable software that can accommodate a large number of atoms interacting through complex forcefields. To address this, we developed and implemented strategies in the GENESIS molecular dynamics package designed for large numbers of processors. Long-range electrostatic interactions were parallelized by minimizing the number of processes involved in communication. A novel algorithm was implemented for nonbonded interactions to increase single instruction multiple data (SIMD) performance, reducing memory usage for ultra large systems. Memory usage for neighbor searches in real-space nonbonded interactions was reduced by approximately 80%, leading to significant speedup. Using experimental data describing physical 3D chromatin interactions, we constructed the first atomistic model of an entire gene locus (GATA4). Taken together, these developments enabled the first billion-atom simulation of an intact biomolecular complex, achieving scaling to 65,000 processes (130,000 processor cores) with 1 ns/day performance. Published 2019. This article is a U.S. Government work and is in the public domain in the USA.

Optimal Temperature Evaluation in Molecular Dynamics Simulations with a Large Time Step.

In molecular dynamics (MD) simulations, an accurate evaluation of temperature is essential for controlling temperature as well as pressure in the isothermal-isobaric conditions. According to the Tolman's equipartition theorem, all motions of all particles should share a single temperature. However, conventional temperature estimation from kinetic energy does not include Hessian terms properly, and thereby, the equipartition theorem is not satisfied with a large time step. In this paper, we show how to evaluate temperature the most accurately without increasing computational cost. We define two kinds of kinetic energies, evaluated at full- and half-time steps that underestimate or overestimate temperature, respectively. A combination of these two kinetic energies provides an optimal instantaneous temperature up to the third order of the time step. The method is tested for a one-dimensional harmonic oscillator, pure water molecules, a Bovine pancreatic trypsin inhibitor (BPTI) protein in water molecules, and a hydrated 1,2-dispalmitoyl- sn-phosphatidylcholine (DPPC) lipid bilayer in water molecules. In all tests, the optimal temperature estimator fulfills the equipartition theorem better than existing methods and reproduces well the usual physical properties for time steps up to and including 5 fs.

Acceleration of cryo-EM Flexible Fitting for Large Biomolecular Systems by Efficient Space Partitioning.

Flexible fitting is a powerful technique to build the 3D structures of biomolecules from cryoelectron microscopy (cryo-EM) density maps. One popular method is a cross-correlation coefficient-based approach, where the molecular dynamics (MD) simulation is carried out with the biasing potential that includes the cross-correlation coefficient between the experimental and simulated density maps. Here, we propose efficient parallelization schemes for the calculation of the cross-correlation coefficient to accelerate flexible fitting. Our schemes are tested for small, medium, and large biomolecules using CPU and hybrid CPU + GPU architectures. The scheme for the atomic decomposition MD is suitable for small proteins such as Ca2+-ATPase with the all-atom Go model, while that for the domain decomposition MD is better for larger systems such as ribosome with the all-atom Go or the all-atom explicit solvent models. Our methods allow flexible fitting for various biomolecules with reasonable computational cost. This approach also connects high-resolution structure refinements with investigation of protein structure-function relationship.

Kinetic energy definition in velocity Verlet integration for accurate pressure evaluation.

In molecular dynamics (MD) simulations, a proper definition of kinetic energy is essential for controlling pressure as well as temperature in the isothermal-isobaric condition. The virial theorem provides an equation that connects the average kinetic energy with the product of particle coordinate and force. In this paper, we show that the theorem is satisfied in MD simulations with a larger time step and holonomic constraints of bonds, only when a proper definition of kinetic energy is used. We provide a novel definition of kinetic energy, which is calculated from velocities at the half-time steps (t - Δt/2 and t + Δt/2) in the velocity Verlet integration method. MD simulations of a 1,2-dispalmitoyl-sn-phosphatidylcholine (DPPC) lipid bilayer and a water box using the kinetic energy definition could reproduce the physical properties in the isothermal-isobaric condition properly. We also develop a multiple time step (MTS) integration scheme with the kinetic energy definition. MD simulations with the MTS integration for the DPPC and water box systems provided the same quantities as the velocity Verlet integration method, even when the thermostat and barostat are updated less frequently.

GENESIS 1.1: A hybrid-parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms.

GENeralized-Ensemble SImulation System (GENESIS) is a software package for molecular dynamics (MD) simulation of biological systems. It is designed to extend limitations in system size and accessible time scale by adopting highly parallelized schemes and enhanced conformational sampling algorithms. In this new version, GENESIS 1.1, new functions and advanced algorithms have been added. The all-atom and coarse-grained potential energy functions used in AMBER and GROMACS packages now become available in addition to CHARMM energy functions. The performance of MD simulations has been greatly improved by further optimization, multiple time-step integration, and hybrid (CPU + GPU) computing. The string method and replica-exchange umbrella sampling with flexible collective variable choice are used for finding the minimum free-energy pathway and obtaining free-energy profiles for conformational changes of a macromolecule. These new features increase the usefulness and power of GENESIS for modeling and simulation in biological research. © 2017 Wiley Periodicals, Inc.

Multiple program/multiple data molecular dynamics method with multiple time step integrator for large biological systems.

Parallelization of molecular dynamics (MD) simulation is essential for investigating conformational dynamics of large biological systems, such as ribosomes, viruses, and multiple proteins in cellular environments. To improve efficiency in the parallel computation, we have to reduce the amount of data transfer between processors by introducing domain decomposition schemes. Also, it is important to optimize the computational balance between real-space non-bonded interactions and reciprocal-space interactions for long-range electrostatic interactions. Here, we introduce a novel parallelization scheme for large-scale MD simulations on massively parallel supercomputers consisting of only CPUs. We make use of a multiple program/multiple data (MPMD) approach for separating the real-space and reciprocal-space computations on different processors. We also utilize the r-RESPA multiple time step integrator on the framework of the MPMD approach in an efficient way: when the reciprocal-space computations are skipped in r-RESPA, processors assigned for them are utilized for half of the real-space computations. The new scheme allows us to use twice as many as processors that are available in the conventional single program approach. The best performances of all-atom MD simulations for 1 million (STMV), 8.5 million (8_STMV), and 28.8 million (27_STMV) atom systems on K computer are 65, 36, and 24 ns/day, respectively. The MPMD scheme can accelerate 23.4, 10.2, and 9.2 ns/day from the maximum performance of single-program approach for STMV, 8_STMV, and 27_STMV systems, respectively, which correspond to 57%, 39%, and 60% speed up. This suggests significant speedups by increasing the number of processors without losing parallel computational efficiency. © 2016 Wiley Periodicals, Inc.

Biomolecular interactions modulate macromolecular structure and dynamics in atomistic model of a bacterial cytoplasm.

Biological macromolecules function in highly crowded cellular environments. The structure and dynamics of proteins and nucleic acids are well characterized in vitro, but in vivo crowding effects remain unclear. Using molecular dynamics simulations of a comprehensive atomistic model cytoplasm we found that protein-protein interactions may destabilize native protein structures, whereas metabolite interactions may induce more compact states due to electrostatic screening. Protein-protein interactions also resulted in significant variations in reduced macromolecular diffusion under crowded conditions, while metabolites exhibited significant two-dimensional surface diffusion and altered protein-ligand binding that may reduce the effective concentration of metabolites and ligands in vivo. Metabolic enzymes showed weak non-specific association in cellular environments attributed to solvation and entropic effects. These effects are expected to have broad implications for the in vivo functioning of biomolecules. This work is a first step towards physically realistic in silico whole-cell models that connect molecular with cellular biology.

Graphics Processing Unit Acceleration and Parallelization of GENESIS for Large-Scale Molecular Dynamics Simulations.

The graphics processing unit (GPU) has become a popular computational platform for molecular dynamics (MD) simulations of biomolecules. A significant speedup in the simulations of small- or medium-size systems using only a few computer nodes with a single or multiple GPUs has been reported. Because of GPU memory limitation and slow communication between GPUs on different computer nodes, it is not straightforward to accelerate MD simulations of large biological systems that contain a few million or more atoms on massively parallel supercomputers with GPUs. In this study, we develop a new scheme in our MD software, GENESIS, to reduce the total computational time on such computers. Computationally intensive real-space nonbonded interactions are computed mainly on GPUs in the scheme, while less intensive bonded interactions and communication-intensive reciprocal-space interactions are performed on CPUs. On the basis of the midpoint cell method as a domain decomposition scheme, we invent the single particle interaction list for reducing the GPU memory usage. Since total computational time is limited by the reciprocal-space computation, we utilize the RESPA multiple time-step integration and reduce the CPU resting time by assigning a subset of nonbonded interactions on CPUs as well as on GPUs when the reciprocal-space computation is skipped. We validated our GPU implementations in GENESIS on BPTI and a membrane protein, porin, by MD simulations and an alanine-tripeptide by REMD simulations. Benchmark calculations on TSUBAME supercomputer showed that an MD simulation of a million atoms system was scalable up to 256 computer nodes with GPUs.

Dimensionality of Collective Variables for Describing Conformational Changes of a Multi-Domain Protein.

Collective variables (CVs) are often used in molecular dynamics simulations based on enhanced sampling algorithms to investigate large conformational changes of a protein. The choice of CVs in these simulations is essential because it affects simulation results and impacts the free-energy profile, the minimum free-energy pathway (MFEP), and the transition-state structure. Here we examine how many CVs are required to capture the correct transition-state structure during the open-to-close motion of adenylate kinase using a coarse-grained model in the mean forces string method to search the MFEP. Various numbers of large amplitude principal components are tested as CVs in the simulations. The incorporation of local coordinates into CVs, which is possible in higher dimensional CV spaces, is important for capturing a reliable MFEP. The Bayesian measure proposed by Best and Hummer is sensitive to the choice of CVs, showing sharp peaks when the transition-state structure is captured. We thus evaluate the required number of CVs needed in enhanced sampling simulations for describing protein conformational changes.